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Abstract 


This  work  examines  the  performance  of  the  Recursive  Least  Squares  (RLS)  Algorithm 
and  a  derivative,  the  Recursive  Lattice  Least  Squares  (RLLS)  Algorithm  in  matching  a 
linear  model  to  a  simple  nonlinear  system.  The  response  of  a  single  degree  of  freedom 
mass-spring-dashpot  system  to  continuous  forcing  is  simulated,  and  estimates  for  the 
modal  parameters  are  obtained.  Nonlinearity  is  introduced  by  allowing  the  restoring 
spring  to  slide  without  friction  in  a  gap  of  specified  width.  Such  a  nonlinearity  is  of 
interest  since  it  is  a  simple  model  of  the  effect  of  loose  joints  in  a  deployable  spacecraft 
structure.  The  effect  of  the  gap  size  on  the  algorithms  is  studied,  and  the  noise  level  and 
sampling  rate  are  adjusted  to  see  how  they  influence  these  results.  The  RLS  Algorithm 
is  found  to  be  the  most  reliable.  Finally,  the  RLS  algorithm  is  extended  to  include  an 
estimate  for  gap  width,  and  good  results  are  obtained  for  the  noise-free  case.  Additional 
conclusions  are  made  on  how  to  minimize  the  effect  of  such  a  gap  nonlinearity  on  the 
identification  process. 
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Notation 


System  Notation: 

A  average  displacement  amplitude 

c  damping  constant 

k  spring  stiffness 

M  mass 

u  input  (force) 

y  output  (displacement) 

yo  location  of  boundary  within  gap 

y(0)  initial  displacement 

Modal  Parameters: 

w  natural  frequency 

f  damping  ratio 

M  mass 

8  gap  width 

Algorithm  Variables: 

а,  a*  coefficients  of  outputs  in  ARMAX  form  (RLLS) 

б, 6*  coefficients  of  inputs  in  ARMAX  form  (RLLS) 
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e  errors  in  equations  (RLS) 

e  forward  eqn  errors,  used  as  regressors  (RLLS) 


G 


forward  coefficient  matrix  < 


6r 


(RLLS) 


(RLLS) 

H  unknown  coefficients  of  parameters  (RLS) 

K,  K*  matrix  of  coefficients  to  regressors  r,  e  (RLLS) 
kj  k*  Kalman  gain  factor  (RLLS) 

n  number  of  steps  in  regression  (RLLS) 

P,  P*  covariance  matrix  (RLLS) 

p  desired  parameters  (RLS) 

g-1  backwards  time  shift  (RLLS) 

r  backward  eqn  errors,  used  as  regressors  (RLLS) 

z  known  quantities  (RLS) 

7  parameter  for  variance  estimate  (RLLS) 

A  forgetting  factor  (1.0  for  no  forgetting)  (RLLS) 

(Note  that  starred  variables  indicate  use  in  the  backwards 
difference  equation  in  the  RLLS  algorithm) 


H 


backward  coefficient  matrix 


o*T 

b'T 


Method  of  Averaging  Variables: 

A  The  constant  amplitude  chosen  for  study 

A(t)  displacement  amplitude  as  a  function  of  time 
a,-  small  angle  variables  for  use  in  computing  boundaries 

(3  change  in  phase  angle  as  a  function  of  time 


4>{t)  frequency  as  a  function  of  time 


Chapter  1 


Introduction 


1.1  Background 


As  man’s  experience  with  space  has  grown,  the  objects  placed  into  orbit  have  become 
larger  and  more  complicated.  Current  project  designs  will  require  an  entirely  new  type 
of  structure  to  provide  the  required  surface  area  or  volume.  Large  sensor  arrays,  wide 
mirrors  for  astronomy,  and  truss  structures  like  the  proposed  space  station  all  fall  into 
this  category  of  Large  Space  Structures  (LSS).  Several  NASA  experiments,  including 
the  MAST  flight  experiment  and  the  COFS  (Control  of  Flexible  Structures)  program, 
deal  with  this  this  type  of  structure. 

Accurate  modelling  of  large  spacecraft  structures  presents  many  novel  difficulties.  It 
is  impossible  to  conduct  full  size  ground  tests  of  the  structure  because  it  is  so  large,  and 
any  finite  element  model  is  therefore  of  limited  accuracy.  Perversely,  the  need  for  active 
control  of  the  structured  dynamics  makes  accurate  modelling  necessary.  One  solution 
to  this  problem  is  to  perform  model  or  parameter  identification,  which  when  applied  to 
structural  dynamics  is  known  as  modal  analysis.  This  involves  the  use  of  a  computer 
algorithm  to  transform  the  motion  of  the  structure  into  parameters  that  can  be  used  in  a 
model  for  control  purposes.  Experimental  modal  analysis  of  large  spacecraft  structures 
presents  unique  challenges,  however.  Because  a  structure  in  space  cannot  be  secured, 
it  will  have  several  rigid  body  modes.  A  further  complication  is  that  the  structure 


may  require  the  application  of  control  forces,  even  while  parameter  identification  is  in 
progress.  Lastly,  since  it  is  being  assembled  on  orbit,  it  may  include  relatively  loose 
joint  hardware,  designed  to  enhance  ease  of  assembly.  These  joints  introduce  important 
nonlinearity  into  the  response.  [11] 

If  modal  analysis  is  used  to  analyze  these  spacecraft  structures,  it  must  have  several 
important  characteristics: 

•  Recursive.  The  algorithm  must  be  capable  of  frequently  updating  the  estimates 
of  parameters,  to  allow  the  structure  to  be  modelled  even  in  the  process  of  mod¬ 
ification.  More  importantly,  recursive  algorithms  require  less  data  storage,  since 
they  only  work  with  the  most  current  information.  This  more  than  compensates 
for  then  disadvantage  that  recursive  (or  ‘on-line’)  algorithms  usually  do  not  give 
as  accurate  results  as  batch  (or  ‘off-line)  algorithms. 

•  Allow  Input  Superposition.  To  prevent  rigid  body  motion  or  excessive  deflections 
during  the  tests,  control  forces  would  have  to  be  applied  while  the  test  was  con¬ 
ducted.  Algorithms  that  require  a  specific  forcing  function  (such  as  impulse  or 
sinusoidal)  could  not  be  used. 

•  Robust.  The  algorithm  must  be  able  to  perform  in  the  presence  of  the  nonlinearity 
introduced  by  the  loose  joint  hardware,  and  give  reliable  estimates. 

Recently,  some  time-domain  parameter  identification  algorithms  have  been  proposed 
which  have  the  first  two  characteristics,  but  there  has  been  little  study  of  their  perfor¬ 
mance  in  the  face  of  nonlinearity.  If  these  algorithms  are  to  be  used  with  confidence 
in  a  parameter  identification  process  on  a  large  spacecraft  structure,  this  performance 


must  be  examined. 


1.2  Thesis  Outline 


The  most  common  type  of  nonlinearity  found  in  the  joint  hardware  of  truss  beams  is  a 
gap  (deadband)  nonlinearity.  Therefore,  in  chapter  2,  a  simple  model  containing  a  gap 
nonlinearity  is  chosen  and  analyzed,  and  the  equations  of  motion  are  developed.  This 
nonlinear  system  will  be  used  to  test  the  identification  algorithms. 

Two  Least-Squares  based  algorithms  are  chosen  for  study,  a  basic  Recursive  Least 
Squares  (RLS)  Algorithm,  and  a  Recursive  Lattice  Least  Squares  (RLLS)  Algorithm 
as  formulated  by  Voss.  These  algorithms,  along  with  other  recursive  algorithms,  are 
clearly  explained  in  a  text  by  Ljung  on  recursive  identification.  [9]  The  basic  theory 
of  the  chosen  algorithms  is  explained  in  chapter  3.  Chapter  4  then  gives  details  on 
the  implementation  of  the  model  and  algorithms  as  part  of  the  computer  simulation, 
including  the  choice  of  forcing  function  and  programming  details. 

The  effect  of  the  nonlinearity  on  the  modal  parameter  estimates  (natural  frequency, 
damping  ratio,  and  mass)  are  studied  for  both  algorithms.  The  influence  of  noise, 
sampling  rate,  and  (for  the  RLLS  algorithm)  the  noise  parameter  7  are  also  examined. 
These  are  detailed  in  Chapter  5.  A  linear  relationship  is  found  between  the  size  of  the 
gap  and  the  estimated  parameters.  Noise,  sampling  rate,  and  choice  of  7  are  found  to 
have  a  significant  effect  on  these  estimates.  An  attempt  to  improve  the  identification 
process  is  made  by  extending  the  RLS  algorithm  to  include  an  estimate  for  the  gap, 
and  good  results  are  obtained  for  the  noise-free  case.  The  results  of  this  extension  are 
shown  in  Chapter  6. 

The  RLS  algorithm  is  found  to  be  generally  superior  to  the  RLLS  algorithm.  Knowl¬ 
edge  gained  from  studying  the  effect  of  the  gap  on  the  three  algorithms  leads  to  several 
conclusions  of  how  to  minimize  the  nonlinear  effect  on  the  identification  process. 


Chapter  2 


The  Nonlinear  Model 


2.1  Analytical  Model 


A  truss  structure  assembled  in  space  typically  will  have  a  static  load-deflection  curve 
as  shown  in  Figure  2.1,  indicating  a  gap  (deadband)  nonlinearity  and  slight  hysteresis. 
In  order  to  study  the  effect  of  this  gap  (without  the  hysteresis)  on  recursive  parameter 
identification  (PID)  algorithms,  a  single  degree  of  freedom  mass-spring-dashpot  model 
is  chosen  (figure  2.2). 


The  boundary  of  the  system  is  allowed  to  move  in  a  frictionless  gap  of  arbitrary 
width,  26.  This  yields  the  following  equations  of  motion: 


My  +  cy  +  k{y  -  6)  =  u  (2.1) 

My  -f  cy  +  k(y  +  5)  =  u  (2.2) 

My  =  u  (2.3) 

c(y  -  yo)  +  *(y  -  yo)  =  o  (2.4) 


where  (2.1)  is  for  motion  on  the  right  side  of  the  gap,  (2.2)  is  for  motion  on  the  left 
side  of  the  gap,  and  (2.3-2. 4)  is  for  motion  inside  the  gap.  The  boundaries  between  the 
states  occur  when 

c(y  -  yo)  +  *(y  -  yo)  =  o 
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DEFLECTION 


Figure  2.1:  Typical  Load-Displacement  Curve  for  Truss  Structure. 


Figure  2.2:  Mass-Spring-Dashpot  Model  for  System  with  Gap  Nonlinearity. 


where  the  stop  location  yo  —  6  at  the  right  boundary  and  yo  =  -6  at  the  left  boundary. 
The  sign  of  the  yo  term  indicates  the  direction  the  system  is  moving  across  the  boundary 
condition.  This  model  has  desired  the  type  of  Load-Displacement  curve. 


Since  the  parameters  of  interest  are  the  modal  parameters  of  natural  frequency  (u>), 
damping  ratio  (f),  and  modal  mass  (m),  the  equations  of  motion  are  rewritten  using 
the  identities  c  =  2fo> m  and  k  =  mw2.  This  produces  equations  of  the  form 


y  +  2fwy  +  w2(y  -  6)  =  ti/m 

(2.5) 

y  +  2fwy  +  co2(y  +  6)  =  u/m 

(2.6) 

a 
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(2.7) 

2fw(y  -  y0)  +  wa(y  -  yo)  =  0. 

(2.8) 

2.2  Adaptation  for  Simulation 


In  order  to  make  these  equations  useful  for  computing  a  displacement  history,  they 
are  changed  to  finite  difference  equations  using  the  Central  Difference  Method.  This 
method  is  chosen  because  of  its  simplicity,  stability,  and  lack  of  numerical  damping. 
This  explicit  integration  method  uses  the  following  approximations  for  the  derivatives: 


y  *  yt+At  -  vt-At) 

V  «  (vt+At  -  2 yt  +  yt-At) 


i 


H 

t**i;*T 


So  the  equations  of  motion  may  be  written  explicitly  to  solve  for  y*+A<  in  e&ch  of  the 
cases.  This  gives  the  final  solution  for  the  motion  which  the  computer  can  solve: 


KM 

if 

$ 

»««;» 

r»yi 


[2  -  (u>A*)2]  y*  +  [-1  +  fwAt]  y<_A»  +  u*  +  (wA*)2  £ 
Vt+At  “  (1  +  fwAt) 


[2  -  (wAt)*l  yt  +  [-1  +  ?wAt]  yt-At  +  vu<~  (wAt)1 5  <rtX 

V‘+A*  =  - (T+7 ^ - -  (210) 


V«+A<  =  2  y*  y«-A<  +  u< 

m 


(2.11) 


i 

q 

1 

$ 


(u>At)2 .  \  /  .  *  \ 

yot^A<  =  ^-(y»  -  yot)  +  (yt  +  At  -  yt_A<  +  yo«-At) 


(2.12) 


Note  that  the  first  two  equations  differ  from  the  linear  system  equation  by  the  term 

(wAt)1  S 
(1  +  fwA  t) 


$ 

$ 

ini 

M 

vW 

is 


which  may  be  thought  of  as  a  constant  noise  term,  or  as  a  disturbance  force.  The 
difference  equations  are  used  to  calculate  the  displacement  of  the  system  at  each  time 
step  yt+Ati  based  on  the  forcing  and  displacement  history.  This  generates  the  pseudo¬ 
data  which  is  used  to  test  the  parameter  identification  schemes. 


In  the  actual  algorithm,  a  slightly  different  formulation  of  the  last  equation  was  used. 
Letting  the  spring  stretch  q  =  y-yo,  q  =  y~yo,  where  yo(0)  =  -6  and  y(0)  =  —S— 
for  the  mass  moving  right  in  the  gap,  then  g(0)  =  — 

Since  the  equation  of  motion  is 

c(y  -  y0)  +  *(y  -  yo)  =  0 


then  cq  +  kq  =  0. 


16 


This  yields  q{tt)  =  ?(0)e  (fc/c) **>  or 


»(«,)  = 

This  sh$ws  directly  the  exponential  nature  of  the  relaxation  of  the  spring/dashpot 
element  within  the  gap.  It  does  however  make  the  assumption  that  the  system  returns 
to  equilibrium  before  touchdown  on  the  other  side.  In  retrospect,  working  with  difference 
equations  directly  might  be  preferable.  The  displacement  is  updated  in  the  gap  as 


At2 

yt+At  =  2 yt  -  yt-A«  + 


The  boundary  conditions  occur  when  y  =  yo  -  j(y  —  Vo)-  Liftoff  occurs  from  the  left 
side  when  the  mass  was  moving  to  the  right,  and  yo  =  —8,  yo  =  0: 

yt+M  >-6-  J^(y<  -  Vt-At)- 

Liftoff  occurs  from  the  right  side  when  the  mass  was  moving  to  the  left,  and  yo  =  -H>, 
Sto  =  0: 

2f 

yt+At  <  8  -  -  yt-At)- 

Initial  conditions  on  the  displacement  are  calculated  by  setting  the  initial  displace¬ 
ment  to  y(0)  and  the  initial  force  and  velocity  to  zero.  This  yields: 

^(v(l)  -  Sv(0)  +  v(-i))  +  ^(v(o)  +  8)  =  ^> 

For  y(-l)  ft*  y(l)  (at  a  minimum  point), 

s(l)  =  \ IP  -  „(0)  +  (^  -  "’«)  A'1]- 

The  actual  algorithm  is  included  in  Appendix  C. 
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Chapter  3 


Identification  Algorithms 


3.1  Choice  of  Domain  and  Algorithms 


Identification  algorithms  are  designed  to  take  the  frequency  response  or  time  history 
of  a  system  and  from  this  information  extract  the  parameters  to  fit  a  specific  model. 
A  flowchart  of  this  process  is  shown  in  figure  3.1.  These  parameters  can  then  be  used 
to  verify  or  update  a  model,  improving  the  accuracy,  to  determine  the  effects  of  a 
structural  modification,  or  to  perform  modal  synthesis  (determining  the  response  of  a 
structure  based  on  the  response  of  its  components).  Identification  algorithms  fall  into 
two  categories:  frequency  domain  methods  and  time  domain  methods. 

In  the  frequency  domain,  the  equations  of  motion  for  a  linear  structural  system  are 
typically  given  as 

\-u/2M  +  jurC  +  K]  y{jv)  =  u(ju>) 

or  y(/w)  =  where  H (jw)  is  the  frequency  response  function  (FRF)  matrix. 

Fullekrug  [4]  classifies  the  frequency  domain  algorithms  into  two  categories:  Sine-Dwell 
methods,  which  measure  the  resonance  of  a  structure  directly  at  its  natural  frequencies, 
and  FRF  techniques  which  measure  responses  for  closely  spaced  discrete  frequencies 
and  use  matrix  techniques  to  extract  the  parameters. 


In  the  time  domain,  the  equations  of  motion  may  be  written  as 

My  +  Cy  +  Ky  =  u(t). 

Fullekrug  divides  the  time  domain  algorithms  into  those  based  on  free  vibration,  and 
those  using  forced  vibration. 

Because  for  this  project  a  recursive  algorithm  allowing  force  superposition  was  de¬ 
sired,  the  algorithms  was  chosen  from  the  time  domain  algorithms  based  on  forced 
vibration.  These  include  such  Least-Squares  based  algorithms  as  the  Eigensystem  Re¬ 
alization  Algorithm  (ERA)  and  Autoregressive  Moving  Average  (ARMA)  methods. 

Two  algorithms  were  chosen  for  study,  a  simple  Recursive  Least  Squares  (RLS) 
algorithm,  and  a  Recursive  Lattice  Least  Squares  (RLLS)  algorithm.  The  latter  was 
chosen  primarily  because  of  its  use  for  identification  in  a  concurrent  experiment  by 
Blackwood.  [2] 


3.2  The  RLS  Algorithm 


The  RLS  algorithm  is  simple  in  concept,  since  it  is  just  a  slightly  modified  form  of  the 
basic  Least  Squares  formulation.  The  idea  behind  a  least  squares  solution  is  simple. 
Given  a  set  of  equations  which  linearly  involve  a  set  of  parameters,  one  may  write  them 
in  matrix  form: 


where 

H  =  Results 

=  Coefficient  vector 
9  =  parameter  vector 
et-  =  error. 

If  the  number  of  equations  exceeds  the  number  of  parameters  (i.e.  n  >  m),  the 
estimated  values  for  the  parameters  §  that  produce  the  minimum  sum  of  the  squared 
errors  e ?  may  be  found  from  the  equation 

9  =  (<$T$)_1$rz. 

This  is  referred  to  as  the  pseudoinverse.  The  relation  of  the  pseudoinverse  form  to  a 
kalman  filter  formulation  is  shown  in  Appendix  B. 
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Since  this  algorithm  is  trying  to  fit  the  modal  parameters  of  a  linear  system, 
y  +  2f<uy  +  w2y  =  u/M,  the  governing  equation  is  put  in  the  form 


and  evaluated  at  a  sequence  of  instants  to  generate  a  matrix  equation.  Given  the  pseudo¬ 
data  points  of  ti{  and  y*,  the  values  of  fw,  w2,  and  l/M  are  found  that  minimize  the 
squared  error.  The  derivatives  of  y  are  calculated  using  the  central  difference  method. 
The  algorithm  is  made  recursive  by  evaluating  the  contrib  tions  to  the  3x3  and 
3x1  $Tz  matrices  at  each  time  step  and  adding  them  to  the  existing  values.  The  new 
matrices  are  used  to  calculate  new  values  of  9_  for  each  t  >  3. 

3.3  The  RLLS  Algorithm 

The  RLLS  algorithm  is  a  linear,  discrete  time,  state-space  method.  It  is  based  on 
the  ARMAX  (Autoregressive  Moving  Average,  Exogeneous  inputs)  form  of  the  linear 
difference  equation,  which  may  be  written  as 

y(t)  =  a.\y(t  -  1)  +  6iu(t  -  1)  + - 1-  o„y(t  -  n)  +  bnu(t  -  n)  +  b0u{t)  +  en[t) 

This  form  is  also  known  as  an  estimator  form,  since  the  output  is  a  linear  combination 
of  previous  inputs  and  outputs,  and  may  be  used  to  estimate  future  output  values  for  a 
system,  as  was  done  in  simulating  the  nonlinear  system  (see  equations  2.9  to  2.12).  The 
ARMAX  form  has  been  used  for  many  years  on  linear  systems  with  excellent  results, 
and  the  technique  is  commonly  known  as  Linear  Regression. 

The  RLLS  algorithm  offers  several  advantages.  Besides  being  recursive  and  allowing 
input  superposition,  it  has  fast  convergence  for  small  data  variance,  and  it  can  handle 


many  degrees  of  freedom  simultaneously.  A  disadvantage  is  that  it  requires  strong 
excitation  of  the  natural  frequencies,  and  so  has  difficulty  with  a  pseudorandom  input. 


The  RLLS  algorithm  is  unusual  because  it  also  makes  use  of  a  backwards  difference 
ARMAX  equation: 


y(t  -  n)  =  a^y[t)  +  f>S“{0  + - h  lV(*  “  »  +  1) 

+6n-iu(*  -  n  +  1)  +  6*u(t  -  n)  +  rn(t). 


The  two  equations  then  go  through  a  change  of  basis  to  make  the  error  terms  the 
regressing  coefficients.  This  derivation  is  included  in  Appendix  A.  The  result  is  a  pair 
of  equations  which  form  a  lattice  (see  figure  3.2): 


'  »(<) ' 

l  «w , 

/  \ 

y (t  -  n) 

^  u(t  -  n)  j 


( 


&)(*  ~  1) 
El(*-  1) 


K^t-n+1)  K'n_x{t) 


V  *)  j 

§o(,t  -  n  +  1) 
e1(t-n  +  2) 


V  6*-l  (0  ) 


where  A,  and  A ■  are  2x2  Reflection  Coefficient  matrices  and  e,-  and  r,  are  2x1  error 
vectors  used  as  regressors. 


The  Voss  formulation  of  the  algorithm  makes  use  of  two  Kalman  filters  to  solve  for 
the  values  of  A,  and  A,*.  This  reduces  the  number  of  computations  required.  According 
to  Ljung  [9],  this  type  of  formulation  produces  assymptotic  convergence  of  the  estimates 
with  a  zero  mean,  normal  distribution,  covariance  matrix  for  the  errors  e  and  r.  He  also 
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mentions  that  while  it  works  well  using  a  sinusoidal  input,  it  may  be  bad  for  white  noise 
input.  The  Voss  RLLS  algorithm  in  algebraic  form  is  included  at  the  end  of  Appendix 
A. 


Once  the  K{  and  A,*  have  been  determined,  the  desired  modal  parameters  may  be 

i  , 

•• 

extracted  by  examining  the  equations  of  motion  in  estimator  form  as  computed  through 
the  Central  Difference  Method: 


2^(yt+A*  -  2y«  +  yt-At)  +  ^-(yt+ a«  -  yt-At)  (f^o)  +  y t  (& 


This  yields  the  estimator  form: 


y<+At  — 


2  -  w^A t2 


1  +  fwoAt 


yt  + 


A  t2/M 


1  +  fwoA  t 


ut  + 


The  parameters  may  then  be  extracted,  since 


-1  +  fwpAt' 
.  1  +  £WqA t  . 


y<-A<- 


y«+At  —  GiVt  +  G%  ut  +  G  syt-At 


fwAt 

Wo 

f 

M 


1  +  C» 

1  -Gs 

—  2  -  Gi(l  +  fwAt) 

fwA  t 
wA  t 
6t2/G2 

1  4-  fwA  t 


The  RLLS  algorithm  has  been  studied  in  several  papers,  mostly  by  Morf  for  linear 
models.  [7.8]  This  algorithm  has  a  number  of  advantages  over  the  simpler  RLS  algorithm. 
Unlik-  iLS  algorithm,  with  which  there  is  a  penalty  in  accuracy  for  going  beyond 
the  true  mode’  >rder,  the  RLLS  algorithm  may  be  based  on  a  model  of  arbitrary  order. 
Additionally,  all  estimates  provided  by  lower  order  models  may  be  recovered,  because 
only  the  new  information  provided  by  each  increase  in  model  order  is  used  to  change  the 
estimates.  The  ability  to  overspecify  the  order  enables  an  unknown  system  to  be  given 
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an  arbitrary  order  (baaed  on  an  overspecified  number  of  degrees  of  freedom)  and  a  good 
model  obtained  by  seeing  when  the  results  cease  to  improve  for  a  higher  order.  Another 
advantage  for  selecting  the  RLLS  algorithm  is  that  it  it  more  computationally  efficient 
for  large  models  (n>5),  since  the  computing  cycle  time  is  linear  with  the  number  of 
states,  as  opposed  ,to  a  quadratic  relationship  for  the  RLS.  Since  our  model  system  is 
small,  this  advantage  will  not  be  shown  here. 

A  final  incentive  for  choosing  this  algorithm  for  study  was  its  use  for  parameter 
identification  in  a  concurrent  modal  synthesis  experiment  being  conducted  at  M.I.T.  by 
Blackwood. [2]  This  experiment  involves  the  analysis  of  a  real  vibrating  jointed  beam 
with  a  gap  nonlinearity.  The  bending  modes  of  the  beam  are  excited  by  a  shaker 
apparatus,  and  screws  at  the  joints  may  be  adjusted  to  introduce  varying  amounts  of 
gap  into  the  joint.  The  experimental  beam  is  shown  in  Figure  3.3.  The  spectrum 
analyzer  for  this  experiment,  a  Signology  SP-20,  uses  the  RLLS  algorithm. 
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Chapter  4 


Organization  of  Project 


4.1  Choice  of  Forcing  Function 


One  of  the  most  important  criteria  when  choosing  a  forcing  function  is  that  it  is  ‘per¬ 
sistently  exciting’,  or  excites  all  the  modes  that  need  to  be  identified.  For  his  use  of 
the  RLLS  algorithm,  Morf  generally  assumes  that  input  is  white  noise.  [7]  The  first 
forcing  function  chosen  to  excite  the  system  was  therefore  a  random  input.  This  func¬ 
tion,  however,  did  not  excite  the  system  enough  to  give  good  results,  and  the  parameter 
estimates  given  by  the  parameter  identification  algorithms  were  different  for  each  ran¬ 
dom  function.  Additionally,  Soucy  and  Vigneron  [11]  recommend  against  use  of  random 
forcing  as  tending  to  eliminate  the  nonlinearity  through  averaging. 

Therefore,  the  forcing  function  chosen  for  this  work  is  a  sinusoid  in  a  linear  envelope, 
repeated  six  times  (See  Figure  4.1).  The  basic  function  was  determined  by  Voss  [12] 
to  be  optimal  for  the  RLLS  frequency  estimate,  and  it  is  repeated  to  give  the  mass 
and  damping  ratio  estimates  time  to  converge.  Too  short  a  time  makes  the  algorithms 
provide  meaningless  results  for  larger  gaps,  and  the  RLLS  algorithm  becomes  very 
sensitive  to  choice  of  7.  It  was  determined  by  trial  that  long  term  forcing,  rather  than 
free  decay,  was  needed  to  give  good  estimates.  This  repetition  gives  better  results 
than  a  single  function  with  a  longer  envelope.  This  forcing  function  also  gives  a  fairly 
constant  displacement  amplitude,  which  is  important  when  analyzing  the  influence  of  a 
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Figure  4.1:  Forcing  Function 

gap  nonlinearity.  Ideally,  the  forcing  frequency  is  chosen  at  the  natural  frequency  of  the 
system,  but  is  chosen  here  at  95%  of  this  value  to  simulate  modelling  error,  or  inexact 
response  from  actuators.  Both  the  input(force)  and  output  (displacement)  at  each  time 
step  were  made  available  to  the  identification  algorithms. 


4.2  Implementation  of  Algorithms 

The  RLS  algorithm  is  simple  to  implement.  The  3x3  matrix  is  inverted  analytically, 
allowing  the  contribution  at  each  time  step  to  be  added  directly,  rather  than  having 
to  invert  the  matrix  at  each  time  step.  The  algorithm  will  not  converge  exactly  in  3 
time  steps,  even  for  a  linear  noise-free  system,  because  the  Central  Difference  Method 
includes  terms  for  and  y*_At  hut  does  not  account  for  the  force  at  the  previous  time 
step,  ut_At,  that  relates  the  two  displacements. 

The  RLLS  algorithm  required  considerably  more  time  to  implement.  Major  problems 
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were  caused  by  numerical  instability  in  the  form  of  roundoff  error.  It  was  determined 
that  a  change  in  the  4th  decimal  place  of  the  input  would  change  the  first  decimal 
place  of  the  parameters.  The  problem  was  eased  by  using  a  UDUr  formulation  for 
the  covariance  matrix  P,  which  factors  the  matrix  into  two  triangular  matrices  and  a 
diagonal  matrix  as  suggested  in  references  [9]  and  [12].  This  formulation  reduces  the 
number  of  addition  operations  in  the  matrix  multiplications,  reducing  the  chance  for 
roundoff  error. 

Another  cause  of  delay  was  that  the  best  sizes  of  the  required  values  for  the  param¬ 
eters  P(0)  and  7  were  ambiguously  treated  in  the  literature.  P(0)  is  the  initial  value 
for  the  covariance  matrix,  which  requires  large  values  (see  Appendix  A).  Values  greater 
than  200000  x  yielded  negligible  improvement,  so  this  was  chosen.  The  parameter  7 
is  a  scalar  weighting  factor  used  to  account  for  noise  in  the  data.  More  exactly,  7  =  A/3, 
where  A  is  the  forgetting  factor  and  /?  is  the  least  squares  weighting  factor  described  in 
Ljung’s  book.  [9]  The  specified  initial  value,  70  =  AV,  is  the  forgetting  factor  times  the 
nondimensionalized  variance  of  the  data  noise.  The  forgetting  factor  is  used  to  discount 
old  data  if  the  system  is  time  varying.  For  our  constant  system,  this  parameter  should 
be  set  to  1  for  no  forgetting.  For  a  noise-free  case,  7  could  be  set  to  zero,  but  because 
of  numerical  errors,  a  value  of  .000001  was  found  to  be  suitable.  Voss  claims  that  the 
best  results  for  noisy  data  occur  when  V  is  set  to  100  times  the  true  variance. 

Because  of  the  forcing  function  chosen,  it  was  also  found  necessary  to  modify  the  way 
that  the  mass  parameter  was  calculated,  since  a  contribution  from  the  force  at  time  t— At 
was  identified.  The  force  is  changed  from  iq(l/Af)  to  tq(l jM){x)  +  tq_  At(l/M)(l-z) 
(since  u«  »  tq-At)-  This  splits  the  influencing  force  over  the  two  time  steps.  In  this 

G  _  ^  (1  —  x)  At*/Af 

1  +  fwoAt  *  1  +  fwoAt  * 


case,  since 


Then  the  mass  is 


1 

m 

o 

m 

m 


(1  +  fwoAt)(Gj  +  GU) 

This  yielded  the  correct  mass  for  the  linear  noise-free  case. 


To  non-dimensionalize  the  size  of  the  gap,  it  should  be  divided  by  the.total  path 
length  that  the  displacement  covers.  This  was  found  by  writing  a  FORTRAN  program 
to  keep  track  of  the  peak  displacements.  This  program,  along  with  the  system  estimation 
and  parameter  identification  algorithms,  is  included  in  Appendix  C. 
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The  simulation  was  organized  into  three  steps.  First  the  forcing  function  was  com¬ 
puted  at  discrete  time  steps.  This  was  then  used  as  input  into  the  system  to  be  studied, 
and  a  file  of  the  resulting  displacements  was  created.  To  simulate  a  sampling  effect, 
the  displacements  were  recorded  for  only  every  other  time  step.  Last  of  all,  both  the 
force  and  displacement  histories  were  fed  into  the  parameter  identification  algorithms. 
Two  different  versions  of  the  output  could  be  obtained:  either  showing  the  identified 
parameters  at  each  time  step,  to  check  things  like  speed  of  convergence,  or  giving  all 
parameters  at  only  the  last  time  step,  to  save  computing  time  when  running  through 
a  large  number  of  cases.  These  final  values  (after  40  periods)  are  the  ones  used  in  the 
results  section. 


m 
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Chapter  5 

Results 


5.1  Addition  of  the  Gap  Nonlinearity 

The  specific  system  used  for  the  simulations  has  the  following 
parameters: 

k  =  20  N/m  oj  =  2 /sec 

c  —  0.4  N-s/m  f  =  .02 

M  =  5.0  kg  S  =  specified 

y(0)  =-2m  Period  =  r 

This  is  then  forced  by  the  function 

u  =  15(1  -  ,08t)  cos  (l.9t) 

repeated  every  seven  periods,  for  a  total  of  42  periods  (total  time  132  sec). 

The  standard  sampling  rate  is  30  per  period.  It  is  found  that  at  higher  rates  (such 
as  60  per  period)  the  JRLLS  Algorithm  proves  to  be  very  unstable  with  the  parameter 
7,  even  with  the  UDUr  and  Kalman  Filter  formulations.  The  sampling  rate  chosen 
proved  to  be  much  more  stable,  probably  because  of  smaller  numerical  roundoff  errors. 
It  should  be  mentioned  that  the  RLS  algorithm  does  not  suffer  from  these  numerical 
difficulties,  and  its  accuracy  improved  with  the  increased  sampling  rate,  but  not  enough 

t 

to  warrant  using  different  rates  for  each  algorithm. 


The  influence  of  the  gap  nonlinearity  is  studied  by  use  of  the  nondimensional  ratio  of 
6/  A  (gap  size/average  amplitude).  This  particular  ratio  is  chosen  because  it  determines 
the  fr fiction  of  the  displacement  amplitude  that  passes  through  the  gap.  With  this 
choice,  when  the  gap  size  approaches  zero  or  the  amplitude  becomes  very  large,  the 
results  approach  those  of  the  linear  system.  The  Anal  estimates  for  each  of  the  three 
parameters  are  plotted  versus  6/  A  in  figure  5.1. 


10.  12  14. 


a)  Frequency  Estimates 
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Estimates. 


The  RLS  algorithm  shows  a  linear  decrease  in  the  estimates  for  the  natural  frequency 


and  damping  ratio,  and  a  linear  increase  in  the  mass  estimate  as  5/  A  increases.  The 


RLLS  algorithm  gives  the  same  results  for  frequency,  since  the  forcing  has  been  opti¬ 


mized  for  this  estimate.  The  RLLS  mass  and  damping  estimates  are  highly  dependent 


on  the  choice  of  7  for  small  gaps  [6/  A  <  5%),  but  converge  to  the  RLS  estimates  for 


larger  gaps.  Approximate  linear  equations  for  the  RLS  estimates  are: 


oje  =  w(.997- ,54  5/A) 


fe  =  f  (.980  -  1.88  5/  A) 


Mt  =  Af(1.011+  1.54  5/A) 


These  may  be  used  to  adjust  a  linear  model  to  account  for  a  gap.  They  also  could  be 


used  to  estimate  the  size  of  the  gap,  given  a  variation  between  the  linear  model  and  the 


estimated  parameters,  but  surely  there  are  better  ways  to  do  this. 


At  first  it  seemed  that  the  RLLS  estimates  were  more  consistent  with  increasing 


5/  A  because  the  system  was  being  forced  more  closely  to  its  natural  frequency  (i.e.  the 


natural  frequency  had  dropped  to  95%  of  the  linear  value  because  of  the  gap).  However, 


changing  the  forcing  frequency  does  not  change  the  trend.  Examination  of  the  time 


traces  (see  figure  5.2)  shows  that  the  estimates  for  larger  gaps  converge  more  quickly. 


This  faster  convergence  probably  comes  from  what  the  algorithm  sees  as  increased 


variance  in  the  data  relative  to  7,  the  variance  parameter,  since  it  is  noted  that  smaller 


values  of  7  (for  the  same  gap  size)  also  make  the  estimates  converge  more  quickly  (see 


figure  5.3).  This  will  lead  to  an  investigation  into  the  effect  of  noise  on  the  RLLS 


estimates. 


5.2  An  Analytical  Comparison 


The  effect  of  the  gap  nonlinearity  on  the  natural  frequency  of  the  system  may  be  esti¬ 
mated  analytically  using  the  method  of  averaging.  The  system  is  cast  in  the  form 

v  +  ulv  =  /(y,  y), 

where 

/(V,y)  =  -2fw0y  +  Wo<5 

/(y,  y)  =  -2fw0y-wo5 
/(y.y)  =  <«^y 

for  the  right  side,  left  side,  and  inside  of  the  gap  (stipulating  an  unforced  system). 

A  solution  is  assumed  of  the  form 

y(t)  =  A(t)  cos  4>{t) 

4>(t)  =  wot  +  0{t) 
y(t)  =  -A(t)wosin^(t). 

Then,  in  terms  of  <f>, 

f  (A(t)  cos  <p(t),  -  A(t)ojo  sin  <p{t))  =  2fw*  A(t)  sin  +  w*£ 

=  2fw*  A(t)  sin  4>(t)  --  Jq8 
—  Wq  A(t)cos^(t) 

for  the  right  side,  left  side,  and  inside  of  the  gap. 

Assuming  that  A(t)  and  <p(t)  can  be  treated  as  constants  over.one  period,  the  bound¬ 
ary  conditions  are  calculated  in  terms  of  <f>.  For  (f>[t )  =  0,  y(t)  =  A(t)  the  maximum  right 


displacement.  Using  small  angle  assumptions,  the  following  identities  are  determined: 
cos(jt/2  -  <*i)  pa  ai  cos(jt/2  +  oj)  « -oj 
cos(3ir/2  -  as)  «  -as  cos(3*/2  +  a«)  w  a< 
sin((*/2  -  a*)  «  1  sin(*r/2  +  aj)  «  1 
8in(3jr/2  -  as)  »  1  sin(3jr/2  +  a<)  «  1. 

As  the  phase  goes  form  0  to  2jt,  each  liftoff  and  landing  condition  is  encountered. 
Since  the  liftoff  and  landing  moving  left  occur  at  6  —  and  —6: 


and 


2? 


A(i)cos^(t)  s=  y(t)  —  8 - [-A(t)wosin$(t)] 

Wq 


^(0  =  jt/2 


8  +  2f  A(t) 
A(t)  ■ 


A(t)  cos  <f>(t )  =  y(t)  =  -8 


4>{t)  =  v/2  + 


m 


Similarly,  the  liftoff  and  landing  moving  right  occur  at  —8  -  ^ y  and  8 


and 


2?, 


A(t)  cos  <f>{t)  =  y(t)  =  -8 - [-A(t)wosin^(f)] 

Wo 


< P(t )  =  3tt/2  H- 


—8  +  2f  A(t) 

m  ■ 


A{t)  cos  <f>{t)  =  y(t)  =  6 


<f>(t)  =  3*/2  + 


A(tY 
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The  values  for  A(t)  and  /?(<)  may  then  be  calculated  from  the  following  equations: 


A(t)  =  5—  /  f(A(t)  cob  <f>(t)t-A(t)uo  ain  4>(t))  *in  j(t)  d<j>  (5.1) 

2ff(do  Jo 

Kl)  =  2n^QA(t)  So  /(AWCOB^(t)»“‘AWMo«n^(<))COB^(0  d<t>  (5-2) 

If  t  is  chosen  at  a  point  when  A(t)  =  A,  a  constant,  this  yields 


«"  “  2^A< 

/  (2f  wjA(t)  sin  <^(£)  +  oj*S)  cos  4>(t)  d<f> 

Jo 

fie/  2+xnr 

+  /  (woA(f)cos^(£))cos^(i)  d<£ 

r3ir/2+  5777+^  „  „ 

+  /  (2fWQj4(£)sin^(£)  -  ojqS)  cos  <f>(t)  d<j> 

J*/*+7fa 

r3x/2+  An- 

+  /  (tL»o^4(t)cos^(£))cos<^(£)  dtp 

■/s’r/2+i^y+Jf 

+  [  (2fWoA(£)sin^(£) +  Wo5)cos^(£)  ). 

■/iW2+zfo 


The  solution  is  then: 


or  (for  j3(0)  =  0) 


wo$ 

JrA 


The  linear  relation  between  natural  frequency  and  gap  size  can  then  be  shown  by  looking 
at  the  phase  equation: 


m 


w05 

wo  t - —t 

7T  A 
.  1 


Thus  fl  =  w0(l  -  .318  j).  This  has  a  slightly  flatter  slope  than  that  obtained  from  the 
algorithm,  but  still  shows  the  linear  relation. 


5.3  The  Effect  of  Noise  on  a  Linear  Model 


If  the  gap  nonlinearity  ia  indeed  seen  as  adding  noise  to  the  system,  a  study  of  the 
effects  of  noise  on  the  linear  model  is  important.  Noisy  displacement  data  is  produced 
by  adding  a  random  perturbation  to  the  sensed  displacement  at  each  time  step.  The 
maximum  amplitude  of  this  noise  is  proportional  to  the  average  displacement  amplitude, 
Ay  =  N A.  The  proportionality  constant  N  is  referred  to  here  as  the  percent  noise.  The 
noisy  data  is  then  passed  to  the  two  parameter  identification  algorithms. 

This  noise  is  found  to  have  a  very  large  impact,  since  the  estimates  of  Least  Squares- 
based  algorithms  have  bias  errors  with  correlated  noise. [10,12]  The  modal  parameter 
estimates  are  plotted  versus  percent  noise  in  figure  5.4.  The  RLS  algorithm  shows 
quadratic  increase  in  frequency  and  mass  estimates,  and  quadratic  decrease  in  the 
damping  ratio  estimate.  The  RLLS  algorithm  shows  quadratic  increase  in  frequency 
and  damping  estimates,  and  an  inverse  relationship  for  the  mass  estimate.  For  esti¬ 
mates  to  be  within  10%  of  the  true  values,  noise  has  to  be  under  5%  for  the  RLS  and 
under  1%  for  the  RLLS  algorithm.  The  mass  estimate  is  the  limiting  factor  in  both 
cases.  It  is  noted  that  the  addition  of  noise  also  makes  the  algorithms  converge  slightly 
faster,  and  makes  the  RLLS  algorithm  less  sensitive  to  7.  This  is  to  be  expected  if  the 
gap  is  treated  as  an  increase  in  noise,  verifying  this  hypothesis. 

The  two  algorithms  yield  the  same  results  when  the  parameter  7  in  the  RLLS  algo¬ 
rithm  approaches  zero  (for  the  noise- free  case).  This  agrees  with  theory,  and  is  true  as 
long  as  the  initial  estimates  are  zero. [3] 

Adjustment  of  the  parameter  7  can  be  made  to  compensate  for  the  increase  in 
frequency  due  to  noise,  with  very  small  values  being  best  in  the  noise-free  case  and 
larger  values  needed  as  the  noise  is  increased.  Setting  7  =  0  does  not  yield  a  convergent 
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estimate  in  the  noise-free  case  because  of  small  numerical  errors,  so  a  value  of  7  = 
.000001  is  chosen.  Voss’  rule  of  thumb  of  100  times  the  true  variance  seems  to  work 
only  for  percent  noise  under  1%,  since  it  was  found  that  for  1%  noise  the  best  value  is 
.0001  and  for  5%  noise  it  is  4.0  .  The  choice  of  a  larger  value  of  7  allows  more  data  to  be 
processed  before  the  estimates  converge,  allowing  the  algorithm  to  compensate  for  the 
noise.  The  parameter,  however,  is  unable  to  compensate  for  the  very  strong  changes  in 
damping  ratio  and  mass  estimates  with  increasing  noise.  These  estimates  would  have  to 
be  normalized  somehow  for  the  parameter  to  affect  all  the  estimates.  Without  a  priori 
knowledge  of  the  modal  parameters,  this  would  be  difficult. 


*  »ose 


a)  Frequency  Estimates 


b)  Damping  Ratio  Estimates 


o 


c)  Mass  Estimates 


Figure  5.4:  Effect  of  Noise  on  the  Frequency,  Damping  Ratio,  and  Mass  Estimates. 
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5.4  Combined  Effects  of  Noise  and  Nonlinearity 


The  prior  sections  summarize  the  effect  of  measurement  noise  and  nonlinearity  when 
each  is  present.  The  effects  on  the  parameter  estimates  are,  in  part,  in  opposite  direc¬ 
tions.  The  next  step  is  to  determine  the-offect  on  the  estimates  when  both  noise  and 
nonlinearity  are  present.  It  turns  out  that  when  noise  is  added  to  the  nonlinear  data, 
the  curves  for  the  variation  of  parameter  estimates  with  gap  are  shifted  accordingly, 
as  shown  in  figure  5.5.  Here,  only  the  optimal  value  of  *y  for  each  noise  level  is  used 
(.00001  for  linear,  .0001  at  1%,  4.0  at  5%).  Low  noise  levels  produce  the  same  linear 
increase  or  decrease  (higher  noise  levels  have  a  greater  deviation  in  shape,  probably 
due  to  incomplete  convergence).  This  shows  that  the  effects  due  to  the  gap  and  the 
noise  can  be  superimposed.  A  small  amount  of  noise  can  completely  offset  the  drop  in 
frequency  due  to  the  gap. 


Chapter  6 

Extension  of  the  RLS  Algorithm  for  an  Explicit 

Nonlinear  Model 


The  RLS  algorithm  may  be  extended  to  include  an  estimate  for  the  size  of  the  gap  itself. 
This  type  of  extension  will  work  for  any  linearity  that  adds  extra  terms  to  the  equation 

of  motion,  such  as  a  cubic  or  quadratic  nonlinearity.  The  new  matrix  formulation  is 

/ 

..  r  l 

y=  2y  y  u  ±1 

1  1/M 

6(jj 2 

This  is  solved  as  the  RLS  algorithm  was,  through  use  of  the  Central  Difference  Method 
to  compute  the  derivatives,  and  extracting  the  parameter  estimates  through  the  pseu¬ 


doinverse.  The  actual  program  used  is  included  in  Appendix  C.  As  indicated  in  the 

A  A 

equation  above,  it  ignores  the  data  within  the  region  of  the  gap,  -8  <  y  <  8. 


This  algorithm  produces  very  good  results  for  the  noise-free  case,  regardless  of  the 
size  of  the  gap.  In  fact,  the  parameter  estimates  become  slightly  more  accurate  as  the 
gap  increases.  In  the  presence  of  noise,  however,  (since  this  algorithm  is  also  based  on 
a  Least  Squares  formulation)  a  strong  bias  is  introduced.  These  results  are  shown  in 
figure  6.1. 
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Figure  6.1:  Extended  RLS  Estimates  for  Frequency,  Damping  Ratio,  Mass,  and  Gap 
for  Three  Gap  Sizes.  51 


Chapter  7 


Conclusions  and  Recommended  Future  Work 

7.1  Conclusions 

A  gap  nonlinearity,  such  as  would  be  found  in  an  assembled  truss  structure,  has  a  very 
definite  effect  on  parameter  identification  algorithms.  Increasing  the  gap  size  produces 
a  linear  decrease  in  frequency  and  damping  ratio  estimates,  and  a  linear  increase  in  the 
mass  estimate  for  both  the  RLS  and  RLLS  algorithms.  The  gap  does  not  cause  the 
algorithms  to  fail,  because  both  simply  fit  the  best  linear  model  they  can.  Additionally, 
the  gap  causes  faster  convergence  of  the  estimates,  and  makes  the  RLLS  algorithm  less 
sensitive  to  choice  of  the  noise  variance  parameter  7. 

The  addition  of  noise  to  the  data  also  produces  faster  convergence  and  less  sensitivity 
to  7.  Noise  however  causes  very  large  bias  errors  in  the  estimates,  and  these  are  in  the 
opposite  direction  to  the  change  caused  by  a  gap.  These  large  biases  are  probably  not 
seen  in  the  results  for  the  gap  nonlinearity  (figure  5.1)  because  the  gap  is  seen  only  as 
a  low  frequency,  constant  amplitude  disturbance.  While  the  Central  Difference  Method 
used  in  the  algorithms  differentiates  the  noise,  making  it  worse,  the  constant  amplitude 
error  from  the  gap  differentiates  to  zero.  Since  bias  errors  for  the  gap  are  probably  small, 
the  drop  in  frequency  and  damping  ratio  noted  with  increasing  gap  should  accurately 
reflect  the  average  system  parameters.  Use  of  the  method  of  averaging  showed  this  same 
type  of  linear  decrease  in  frequency. 


The  effect  of  5%  noise  on  the  three  algorithms  is  shown  in  Table  7.1.  The  optimal 
values  of  the  variance  parameter  7  are  used  throughout.  The  frequency  estimates  of  the 
ERLS  algorithm  in  the  presence  of  noise  are  consistently  worse  than  those  of  the  RLS 
and  RLLS  algorithms,  but  the  damping  and  mass  estimates  are  better  than  those  from 
the  RLLS  algorithm,  and  become  better  than  those  from  the  RLS  algorithm  with  higher 
gap  size.  The  RLS  and  RLLS  results  are  comparable  in  the  magnitude  of  the  frequency 
estimate  error,  but  the  RLS  yields  far  better  estimates  for  the  other  parameters.  In 
general,  for  moderate  noise  levels,  the  RLS  algorithm  is  the  best  choice  of  these  three. 
For  low  noise  levels  and  significant  gap,  the  ERLS  Algorithm  would  be  preferable. 


Table  7.1:  Comparison  of  the  Effect  of  5%  Noise  on  the  Algorithms.  Given 
as  a  percent  change  from  the  noise-free  value  of  the  estimates. 

An  additional  reason  that  the  RLS  algorithm  is  preferable  to  the  RLLS  algorithm 
is  that  the  RLLS  algorithm  is  too  sensitive  to  the  noise  variance  weighting  factor  and 
to  finite  roundoff  errors.  Therefore,  certain  gaps  sizes  do  not  mix  well  with  choices 
of  7,  as  was  found  for  the  choice  of  7  =  .0001  at  the  1%  noise  level  with  a  4.3% 
gap.  Numbers  around  this  unstable  point  yielded  worthless  answers  as  the  algorithm 
could  not  converge.  While  a  faster  sampling  rate  improved  the  performance  of  the  RLS 


algorithm,  it  only  aggravated  the  problems  in  the  RLLS  algorithm. 

The  noise  problem  can  be  dealt  with  in  a  number  of  ways.  The  easiest  way  is  to 
use  one  or  more  low  pass  filters,  to  eliminate  the  noise  before  it  enters  the  algorithms. 
Another  common  method  is  to  overspecify  the  number  of  degrees  of  freedom  in  the 
RLLS  model  for  the  system,  allowing  the  algorithm  to  fit  the  noise  to  high  frequency 
modes.  In  this  case,  the  RLLS  algorithm  may  work  better,  since  the  system  is  easily 
overspecified  and  the  improvement  with  each  extra  degree  of  freedom  may  be  examined 
without  running  the  algorithm  again.  There  is  an  entire  class  of  Least-Squares  based 
algorithms  that  are  designed  to  handle  the  bias  due  to  noise  by  increasing  the  dimensions 
of  the  matrices  to  include  the  magnitude  of  the  error  due  to  noise.  These  are  referred  to 
as  Extended  Least  Squares  algorithms.  Finally,  it  is  probably  better  to  use  acceleration 
data  and  integrate  for  the  other  variables,  thereby  minimizing  the  noise  effects  rather 
than  increasing  them. 

If  a  detailed  model  including  the  gap  is  desired,  the  ERLS  algorithm  as  written  in 
this  thesis  may  be  used,  or  another  algorithm  may  be  similarly  modified.  This  gives  very 
accurate  results  in  the  absence  of  noise.  There  are  also  a  number  of  algorithms  which 
use  different  approaches  to  get  a  value  for  the  nonlinearity  directly.  The  Kalman  filter 
may  be  extended  to  include  a  nonlinearity.  In  the  Bayesian  form  [9],  the  parameter  is 
thought  of  as  a  random  variable,  and  a  solution  is  given  as  a  probability  density  function 
for  the  states.  Billings  [I]  expands  the  nonlinearity  into  sines  and  cosines  and  obtains 
a  linear  model.  Others,  such  as  the  Threshold  Adaptive  Modelling  Technique  [5]  or 
Volterra  and  Wiener  Series  use  a  complicated  kernal  formulation  and  are  quite  recent. 
They  seem  to  show  promise  when  used  with  complex  nonlinearities  like  a  gap. 

Short  of  using  these  more  complicated  techniques,  there  are  several  ways  thst  the 
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effect  of  the  gap  on  the  identification  process  can  be  minimized.  It  is  important  to 
use  a  continuous  forcing  function,  since  in  free  decay  the  size  of  the  gap  relative  to 
the  average  displacement  amplitude  is  always  growing,  leading  to  constantly  greater 
distortion  of  the  results.  The  frequency  content  should  excite  all  the  desired  modes,  and 
the  amplitude  of  th&'-displacement  should  be  relatively  constant.  While  it  would  be  nice 
to  force  the  displacement  amplitude  to  be  as  large  as  possible  (since  this  minimizes  the 
numerical  problem  and  the  nondimension al  gap  and  noise  sizes),  it  is  preferable  that  the 
displacement  amplitude  matches  that  which  the  system  is  expected  to  encounter  during 
use.  In  this  case,  the  linear  nature  of  the  frequency  drop  may  be  exploited  in  that  only 
two  or  three  amplitudes  need  to  actually  be  tested,  and  a  linear  plot  of  frequency  versus 
displacement  amplitude  can  be  drawn.  This  would  allow  a  switch  between  models  at 
different  amplitudes  for  best  results.  The  forcing  needs  to  continue  for  a  long  enough 
time  (here  42  periods  was  chosen,  and  this  is  probably  a  minimum)  for  the  algorithm 
to  converge. 

It  may  also  be  noted  that  the  higher  the  damping  ratio,  the  more  accurate  the 
damping  ratio  estimate.  This  was  found  in  a  separate  project,  not  discussed  here.  The 
effects  of  a  gap  can  be  eliminated  entirely  by  loading  the  structure  in  compression  (or 
tension),  and  never  allowing  the  displacement  to  cross  the  region  of  the  gap.  This  is 
sometimes  done  in  practice. 

With  precautions,  parameter  identification  of  a  system  with  a  gap  nonlinearity  can 


be  conducted,  and  very  good  results  obtained. 


7.2  Recommended  future  work 


This  study  of  the  effect  of  a  gap  nonlinearity  on  parameter  identification  could  be 
followed  up  in  a  number  of  ways. 

1.  Higher  Degree  of  Freedom  Model.  A  model  consisting  of  two  or  more  degrees  of 
freedom,  each  containing  a  gap  nonlinearity  could  be  studied.  This  work  could 
be  related  to  the  modal  synthesis  of  the  Blackwood  experiment.  It  could  be 
determined  whether  it  is  better  to  include  the  nonlinearity  as  a  disturbance  force 
J  —  kS  instead  of  a  change  in  displacement.  Additionally,  a  better  sense  of  how 
the  gaps  interact  over  the  entire  structure  would  be  obtained. 


2.  Reliability  With  Other  Nonlinearities.  A  model  containing  hysteresis,  as  well  as  a 
gap,  or  a  model  with  a  cubic  nonlinearity  included  could  be  tested.  This  should 
give  a  better  indication  of  how  the  algorithms  would  work  in  actual  use. 

3.  Use  on  Real  Structure.  The  algorithms  need  to  be  tested  on  a  real  nonlinear  struc¬ 
ture,  to  determine  the  accuracy  of  the  estimates  compared  to  those  obtained  by 
other  test  methods. 

4.  Study  of  Noise.  Some  of  the  ways  to  minimize  the  effect  of  noise  on  the  algorithms 
could  be  studied,  including  the  addition  of  noise  modes  or  the  use  of  an  Extended 
Least  Squares  formulation  of  the  algorithms.  Further  study  on  the  effect  of  differ¬ 
ent  types  of  noise  could  be  done  analytically,  as  Longman  has  done  for  the  ERA 


algorithm. [10]  The  interaction  of  a  Minimum  Model  Error  (MME)  state  estimator 
[7]  with  noisy  gap  data  would  be  another  interesting  line  of  work. 

5.  Extension  of  Studied  Algorithms.  The  RLLS  algorithm  could  be  extended  to  in¬ 
clude  an  estimate  for  the  gap,  as  was  done  in  this  thesis  for  the  RLS  algorithm. 
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The  Extended  Least  Squares  form  mentioned  above  could  be  used. 

6.  Study  of  Additional  Algorithms.  Study  of  the  effectiveness  on  a  gap  nonlinear 
system  of  the  ERA,  or  a  nonlinear  identifying  algorithm  could  be  undertaken. 
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Appendix  A 


Derivation  of  the  Voss  RLLS  Algorithm 


A  RLLS  algorithm  is  based  on  the  interaction  of  two  linear  difference  equations  in 
ARMAX  (Auto-Regressive,  Moving  Average,  Exogencous  inputs)  form.  The  first  is 
called  the  forward  difference  equation: 

y(0  =  aiy(i  -  1)  +  &iu(t  -  1)  + - f  a„y (t  -  n)  +  bnu(t  -  n)  +  b0u[t)  +  e„(t).  (A.l) 

The  terms  y(t)  and  u(t)  indicate  outputs  and  inputs  at  time  t  respectively,  so  that  the 
current  output  depends  on  a  linear  combination  of  the  previous  inputs  and  outputs.  The 
variable  n  indicates  the  finite  number  of  previous  time  steps,  and  is  known  as  the  order 
of  regression.  The  term  e„(t)  is  the  error  between  the  linear  approximation  of  order  n 
and  the  actual  y(t).  This  form  of  the  equation  is  also  known  as  a  linear  regression. 

The  second  equation  is  known  as  the  backward  difference  equation: 

y(t-n)  =  aoy(t)  +  6ou(t)  +  •  •  •  +  o*_!y(t  -  n  +  1)  +  —  n  +  1) 

+ b*nu(t  -  n)  +  r„(t). 

The  equations  are  shown  graphically  in  figure  A.l  below.  Since  the  format  of  each 
equation  is  similar,  only  the  forward  difference  equation  is  treated  for  now. 

The  ARMAX  equation  may  be  written  in  a  matrix  form: 

y(*)  =  9T  4>  +  e„(t)  (A. 2) 
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Figure  A.l:  A  Graphical  Interpretation  of  the  RLLS  Difference  Equations. 


y(t  -  1) 
u(*  -  1) 

»(0  ~  <ii  bi  •••  an  bn  bo  "b  en(0 

y(f  -  n) 
u(t  -  n) 


where  9  is  then  the  vector  of  unknown  coefficients  a  and  6  and  $  is  a  vector  of  the 
previous  outputs  and  inputs.  One  way  to  solve  this  type  of  equation  is  through  use 
of  a  Kalman  filter,  which  computes  the  unknown  6  based  on  the  following  state  space 
representation: 

0(t  +  l)  =  9(t) 

y  (0  =  pT(t)  *(*)  +  «(*)• 

This  shows  explicitly  that  9(t)  is  made  up  of  values  which  are  constant  for  all  time  t. 
Note  that  0T(t)‘p(t)  =  <pT(t)$(t)  since  both  variables  are  vectors. 


A  more  general  state  space  equation  is  written  as 


z(t  +  1)  =  F  x(t)  +  G  u(t)  +  w(t) 
y(t)  =  H  z(t)  +  I  u(t)  +  e(t) 


Here,  w[t)  indicates  changes  in  the  states,  and  e(t)  indicates  errors  in  the  output. 


Thia  general  state-space  equation  yields  three  iterative  equations  to  solve  for  *(t)  = 
E( x(t)|y(0),u(0), . , . ,  y(t-n),  u(t-n) )  (an  equation  indicating  that  i(t)  is  an  estimation 
based  on  the  previous  data), 


£(t  +  i)  =  *’(«)*(*)  + £(0  (y(0  -  ff(0*(01 

F(t)P(t)HT(t) 

a  ;  H{t)P(t)HT{t)  +  r2(t) 

p(t  + 1)  =  F(t)P(t)Fr(t)  +  n(t)-M0ff(0^(0ii,T(0- 

This  is  known  as  a  Kalman  filter,  k(t)  is  known  as  the  Kalman  gain.  The  parameters 
r\  and  r2  come  from  the  2x2  matrix  of  variances,  R,  where  =  E(w(t)wT(t)),  r2  = 
E(e(t)eT(t))  and  ru  =  E(w(t)eT(t)).  Since  the  errors  are  independent,  r^2  —  0.  P[t)  is 
then  the  diagonal  covariance  matrix  (the  inverse  of  R).  The  derivation  of  these  equations 
may  be  found  in  Ljung  [9j. 

For  our  specialized  case,  F  =  1,  H  =  <pT[t),  G  =  I  =  0,  x(t)  —  and  the  first 
error  r^t)  =  0.  This  yields  the  following  three  equations: 

*>  +  l)  = 

m  = 

p(t+i)  = 

This  formulation  is  equivalent  to  the  least  squares  solution  if  the  system  is  overdeter- 
mincd,  the  state  is  assumed  constant,  and  no  prior  knowledge  of  the  state  is  assumed. 
Although  this  is  surprising,  given  that  the  least  squares  solution  is  deterministic,  and 
the  Kalrnan  filter  is  a  stochastic  system,  this  is  proved  in  Appendix  B.  The  Kalman 
filter  has  the  advantage  that  the  initial  state  can  be  assigned  a  value,  perhaps  from  a 
model.  This  speeds  convergence. 

The  advantage  of  the  RLLS  formulation  is  that  the  effect  of  changing  the  regression 
order  n  on  the  estimate  is  made  explicit.  This  allows  the  number  of  degrees  of  freedom 


0(0  +  *(0  [y(0  -  <pr(O0(Ol 
my® 

<pT(t)P(t)<p(t)  -I-  r2(t) 
P(t)-k(t)pT(t)P(t) 


in  the  model  to  be  overspecified,  and  the  increase  in  order  examined  to  see  when  no 
further  improvement  occurs  in  the  parameter  estimates.  Additionally,  the  estimates  for 
all  orders  less  than  n  may  be  extracted,  as  opposed  to  the  normal  Kalman  filter,  which 
would  require  recomputing,  all  the  coefficients.  The  resulting  form  has  the  shape  of  a 
ladder  (or  lattice)  as  was  shown  in  Chapter  3. 

In  her  thesis,  Voss  makes  maximum  use  of  the  Lattice  properties  and  of  the  power 
of  the  Kalman  filter  by  using  an  individual  Kalman  filter  for  both  the  forward  and 
backward  equations.  She  also  modifies  the  notation  used,  in  order  to  improve  the 
algorithm’s  performance.  She  expands  the  basic  equation  (A.l)  to  include  the  current 
input  as  well  as  output,  and  defines  a  new  variable  *(t),  where  x(t)  =  [  j fr(t)  u T(t)  ]T. 
This  requires  the  assumption  that  y(t)  is  not  affected  by  u(t),  which  simply  means  that 
there  are  no  instantaneous  effects  on  the  system  displacement.  This  yields  the  matrix 
equation 


(  v(0 
\  «(*) 


y(t  -  1) 
u(t  -  I) 

y(t  -  n) 
u(i  -  n) 

There  is  a  similar  equation  for  the  backward  difference  equation.  This  arrangement  is 
made  so  that  inputs  and  outputs  may  be  paired  at  each  time  step.  No  attempt  is  really 
made  to  estimate  the  force  u(t),  so  the  values  of  c  and  d  are  essentially  meaningless. 
Note  that  although  the  algorithm  is  usable  for  multiple  inputs  and  outputs,  I  will  only 
be  using  it  for  the  single-input,  single-output  case,  so  that  x(t)  =  [  y(t)  u(t)  ]T. 


The  two  RLLS  equations  then  become 


*W(t)  =  <5M(t)  ±M(t  -1)  +  «n(t)  (A.3) 

for  the  forward  direction,  and 

*<">(*  -  n)  =  AW(t)  ^?)(«)  +  r„(t)  (A.4) 

for  the  backward  direction.  The  hat  (')  L  .ates  the  variable  being  estimated.  Note 
that  G  and  H  here  are  not  the  same  as  the  coefficients  in  the  general  state  space 
equation.  Here  the  first  row  of  G  is  6.  The  composition  of  H  is  similar  in  size  and 
meaning.  The  notation  merely  indicates  that  the  recursion  order  is  n,  and  has 

no  effect  on  the  value  of  x(t). 


The  vectors  ^")(t  —  1)  and  are  defined  to  be 


(  *(*  - 1) 

*(0 

±H{t- 1)  = 

*(t  -  2) 

£wW  = 

x(t  -  1) 

^  £(*  ~  n)  j 

,  £(t-n  +  l)  J 

The  superscript  n  continues  to  indicate  the  total  time  steps  needed. 


(A.5) 


The  equations  A.3  and  A.4  then  undergo  a  change  of  basis,  so  that  the  errors  e„(t) 
and  tn(t)  become  the  regressing  terms,  rather  than  z(t).  Since 

6,(0  =  x(n)(t)-G(n)(t)^(t-l)  (A. 6) 

6,(0  =  *<">(*-«) -£(")(*)  £(">(t)  (A. 7) 

so  that 

r^(t  -  1)  =  x W(«  -  „  -  1)  -  £<">(*  -  1  )*W(t  -  1)  (A. 8) 

If  G^(t)  and  are  set  equal  to  a  zero  2x2  matrix  for  all  time  t,  then 


»(<)  =  *,0)«  EoW  =  £|0,(<) 
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which  is  simply  the  data  taken  at  time  t. 


The  equation  A. 8  may  then  be  written  in  matrix  form,  which  Bhows  the  change  of 
basis: 


tod~  !) 

hx2  0 

0 

(  \ 
x(t  -  1) 

“  1) 

= 

-  1)  Iiz2  ••• 

0 

x(t-2) 

4-  ~ff(n-^(t  ~  1)  -» 

hzt 

,  ^  -  n )  , 

=  T^n)(t  -  1) 

which  makes  use  of  the  definition  of  <^n\t  -  1)  in  equation  A.5,  and  also  equation  A.7. 


From  this,  choose  a  matrix  K  partitioned  in  the  form 


K  = 


Kq  Kx  ...  Kn-X 


where  the  Ki  are  the  reflection  coefficient  matrices,  which  are  square  matrices  of  dimen¬ 
sion  equal  to  number  of  rows  in  x  (2  for  a  SISO  system).  Define  K  to  be 

K  —  G  T~l 


so  that  in  the  new  regressors  (from  equation  A. 3), 


xW(0  =  (dW(t)r-1)(r^n)(«-i))  +  &%(0 

/  .  .  \ 


K0  K  i 


Kn-l 


Co(*-  !) 

ri(t  -  l) 


+  «nW 


V  tn-dt  -  1 )  J 

=  £  *»-*(*)  !)+&»(*)• 

This  equation  is  easier  to  work  with,  since  only  one  regressor,  r„(£  —  1),  is  needed 
at  a  time,  making  it  much  more  suitable  for  an  algorithm.  The  covariance  matrix  P 


undergoes  a  similar  change  of  basis  so  that  P  =  TT(t)  P(t)  T{t).  The  new  equation  in 
K  is  solved  by  a  Kalman  filter,  and  becomes  (after  a  shift  in  time) 

l  _  _ ~~  2)rn(*  —  1) _ 

»£(*-  l)fn(<  ~  2)rn(t  -  l)  +  7»(0 
Kn{t)  =  Kn(t-l)+[en{t)- Kn{t-l)rn(t-l)]kT 

Pn(l-l)  =  i[f-*rJ(I-l)]P„(l-2) 

If  the  parameter  identificatior  is  working,  the  values  of  e„(t)  and  will  be  minimized 
as  t  — *  oo  and  these  minimum  values  should  decrease  with  the  order  of  regression  n. 

A  similar  Kalman  filter  is  provided  for  the  backwards  difference  equation,  using  e(t) 
as  the  recursor,  with  all  other  variables  starred.  The  values  for  the  recursors  as  n  is 
increased  are  found  by  several  updating  equations  which  relate  the  regressors  from  both 
sets  of  equations,e(t)  and  r(t)  .  Since 

*(0  =  ]£ffn-l(*)En-l(*-  l)  +  £n(0 

=  X  Kn(t)Ut-  l)  +  &,+1  (0 

VlW  =  «n(0  ~  Knifait  ~  1). 

By  similar  derivation  from  the  backward  Kalman  filter, 

*^+i(0  =  “  1)  ~  K*it)§n(t). 

These  last  two  equations  clearly  show  the  exchange  of  information  between  filters  which 
give  the  lattice  form  shown  in  figure  3.2. 

There  is  also  an  updating  equation  for  the  parameter  7,  which  is  used  to  account 
for  measurement  noise.  The  initial  value  is  7o(t)  =  AV,  the  forgetting  factor  times  the 
assumed  variance  in  the  measurements.  The  updating  equation  is  derived  in  reference 
[9]  from  the  definition  of  the  least  squares  weighting  factor  where  7 n(t)  =  A/?„(t), 

7n+i(0  =  7 »(0  -  A r£(t  -  l)P„(t  -  l)r„(t  -  1). 


Finally,  the  algorithm  calculates  the  unknown  matrix  G  at  each  time  step  and  for 
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each  value  of  n  as  follows.  Given 


i(“+1)(l)  =  £*„(<)r„(«-l) 

=  C<">(l)*<">(l-1)  +  *„(«)  r.(l-l) 


and  that 


x(n+1  )(t)  =  G(n+1)(t)  ^n+1>(t  -  1) 


G^n+1^(<)  ^n+1>(t  -  1)  =  G<n>(i)  ^W(t  -  1)  +  £n(t)rn(t  -  1). 


From  equation  A. 8, 


rjt  -  1)  =  ?<")({  -  n  -  1)  -  £<n)(t  -  1)  £<")(*  -  1) 


so  this  can  be  factored  as 


4>(n\t  -  1) 


G(«+l)W  ^  = 

V  *(*-(»  +  !))  J 

G(rl)(0  4n){t  -  1)  +  Kn(t)  [  -  &(")(<  -  -  1)  +  x(t  -  n  -  1)] 


which  yields 


g("+1)(0  =  (  GW(t)  0  ]  +  Kn{t)  [  -Jf(n)(t  -  1)  I  ]• 


This  simple  updating  equation  requires  much  less  work  than  calculating  G  =  K  T.  The 


updating  equation  for  the  H  matrix  is  similarly  derived  to  be 


H(n+1){t)  =  [  o  tf(")(t  -  1)  ]  +  K{t)  [  /  -GW(i)  ]• 


The  entire  algorithm  is  included  as  an  appendix.  Derivation  of  the  second  Kalman 


filter  for  recursion  in  e„(t)  and  K *  may  be  found  in  Appendix  A  of  [12]. 


An  RLLS  algorithm  then  actually  consists  of  two  loops.  The  inner  loop  estimates 
the  unknown  parameters  as  n  (the  order  of  the  recursion)  increases.  The  outer  loop 
reads  the  next  data  point  as  time  increases  and  updates  in  time  the  two  linear  difference 


Ini tialization: 


X  (0)  ■ 
n 

P  (-1  ) 
n 

V°  ' 


K  (0)  -  0 
n 


P*(0) 

n 


n-O+M-1 
^  - — *  u>o,o  oojr") 


large  I 

XV,  where  X  *  forgetting  factor 

V  m  measurement  noise  variance 


Each  new  measurements 


e  (t)  =»  r  (t)  ■  [ yT  ( t )  u^tt)]1 


For  n=0  +•  M-1 

k  =  — 


P  ( t-2 )  r  (t-1 ) 
n  — n 


r  (t-1 )  P  { t-2 )  r  (t-1 )  +  y  (t) 

— n  n  — n  n 

K  (t)  =  K  (t-1)  +  fe  (t)  -  X  (t-1)  r  (t-1)  Ik1 
n  n  1 — n  n  — n 

P  (t-1)  =  ~  1 1  -  k  rT  (t-1)}P  (t-2) 
n  X  1  —  — n  ‘  n 


* 

k  = 


P  (t-1)  e  (t) 
n  — n 


eT(t)  P*(t-1 )  e  (t)  +  y  (t) 
— n  n  — n  n 


X*(t)  =  X*(t-1)  +  [r  (t-1)  -  X * ( t-1 >  e  (t)]k 
n  n  *•— n  n  — n 

P*(t)  -  1[I  -  k*  eT(t)l  P*(t-D 
n  X 1  —  — n  J  n 


*T 


For  n*M-1 : 

e  , (t)  =  e  (t)  -  X  (t)  r  (t-1) 

— n+1  — n  n  — n 

r  , (t)  =  r  (t-1 )  -  X*(t)  e  (t) 

—n+1  — n  n  — n 

Y  ,(t)  =  Y  (t)  -  X  rT(t-1)  P  (t-1)  r  (t-1) 
n+1  n  — n  n  —n 

To  recover  RLS  coefficients  Oc(t)  «  G(t)  £(t-1)): 

G(n+I)(t)  =  [G(n)(t)  0]  +  X  (t)  [-H(n,(t-D  i] 


H 


(n+1  ) 


(t)  =  [O  H(n5(t-1)]  +  X*(t)  [i  -G(n>(t)] 


with:  G^)(t)  =  ^(t)  and  H<^(t)  =  XQ*(t) 


Figure  A. 2:  The  Voss  Formulation  of  the  RLLS  Algorithm. 
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Appendix  B 

Relationship  Between  the  Pseudoinverse  and 

Kalman  Filter 


This  analysis  is  based  on  a  discussion  in  reference  [3]. 


The  basic  parameter  identification  equation  is 

y(t)  =  <pf  e(t) 

so  that  <pf  and  y(t)  are  the  parameter  coefficient  vector  and  output,  respectively.  Given 
m  time  steps,  the  Kalman  equations  look  as  follows: 


0(t+  1) 

m 

F(t  +  1) 


$(*)  +  *(*)  [HO -*(*)e(01 

$(t)p(t)^r(0  +  Mt) 

P{t)  +  fli(t)  -  k(t)3>r(t)P(t) 


Where 


This  system  has  the  initial  conditions  0(0)  =  0  if  there  are  no  coefficients  determined 
initially.  Since  the  variance  of  the  data  at  t  =  0  is  zero,  the  values  of  the  diagonal 
variance  matrix  are  zero  (i.e.  Fi(0)  =  0  and  /?2(0)  =  0).  At  time  t  —  1,  there  will  be  a 


large  variance  in  the  data  as  the  first  data  point  is  taken,  so  that  R  =  p  x  /  where  p  is 
a  large  number.  Since  the  covariance  matrix  is  P,  where  P  =  i2_1,  P(l)  =  0* 


Starting  at  t  =  0  yields 


e(i) 

m 

o 


£(o)[K(o)] 

P(0)S>T(0) 

$(0)P(C)$T(0) 

P(0)  -  k( 0)  $(0)P(0) 


If  we  set  P(0)  =  P0)  $(0)  =  $0  and  7(0)  =  Vo, 

P0  =  Po^oPoSoV^oPo 
QqPq  =  $oPo$o  ($0Po^o  )~1$oPo 


Which  is  true  for  any  value  of  Po,  so 

0(1)  =  Po^oPo^r'Vo 

Which  is  equivalent  to 

0(1)  =  (<Po^o)_1^oPo7o 

which  is  the  formulation  of  the  problem  in  a  least  squares  solution  with  Po  acting  as  the 
weighting  matrix.  For  this  reason,  Ljung  states  that  this  matrix  is  the  natural  choice 
for  a  weighting  factor.  [9] 

It  is  interesting  that  these  two  very  different  approaches,  one  deterministic  and  the 


other  stochastic,  each  act  to  minimize  the  errors  in  an  analogous  way. 
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Appendix  C 

The  Fortran  Programs 


C.l  The  Nonlinear  System  Model 

PROGRAM  "SYSTEM" 

PROGRAM  MSPR 

COMMON  J . T , XI . X2 . F2 , DT . ROMEGA . RPSI , DISP 
COMMON  RMASS. DELTA, BU.BL.CF 
OPEN (unit=14 ,f ile=f ile , status* ’ old’ ) 

OPEN (unit=15 ,  f ile=f ile , status* ’ old ’ ) 
0PEN(unit*16.file=f ile .status* ’new’) 
0PEN(unit=17.file=f ile .status* ’new’) 

OPEN (unit=50, f ile =f ile .status* ’old’) 

C  INITIALIZATION  OF  PARAMETERS 

romega=2 .0 
rpsi* .02 
rmass=5 .0 

READ(unit=50,fmt=*)  delta 
disp=-2 .0 

READ(unit=14 ,fmt=*)  dt 
cf =2*rpsi/ (romega*dt) 

READ(unit=15,fmt=*)  dummy, f2 
coeff=.5*dt**2 

xl=disp+coef  f  *  (f  2/rmas£»-romega**2*  (disp+delta) ) 
t=0 


x2=disp 


C 

C 


c 

c 


c 

c 


c 

c 


c 

c 


CALCULATING  COEFFICIENTS 
SOT^rpsi^romega+dt 
a=s(2-(romega*dt)  **2)/( 1+SOT) 
b=(SOT-l)/ (1+SOT) 
c=dt**2/rmass/ (1+SOT) 
d= (delta* (romega*dt) *  *2) / (1+SOT) 

BOUNDARY  CONDITIONS 
bu=delta 
bl=-delta 
vel=x2-xl 

IF(vel.GT.O)  GOTO  60 

ON  RIGHT  SIDE  OF  GAP 
20  IF(x2.LT.bu)  GOTO  40 
di8p=a*x2+b*xl+c*f2+d 
CALL  DATA 
GOTO  20 

INSIDE  GAP. MOVING  LEFT 

40  tgap=0 

IF(rpsi.EQ.O)  GOTO  42 

41  z=x2+rpsi*romega*(x2-xl)*exp(-romega/2/rpsi*tgap) 
GOTO  43 

42  z=0 

43  IF (z.LE. -delta)  GOTO  60 
disp=2*x2-xl+c*f 2 

CALL  DATA 
tgap=tgap+dt 
GOTO  41 

ON  LEFT  SIDE  OF  GAP 
60  IF(x2.GT.bl)  GOTO  80 
disp=a*x2+b*xl+c*f2-d 
CALL  DATA 
GOTO  60 


o  o 


\  ' 


INSIDE  GAP, MOVING  RIGHT 

80  tgap=0 

IF(rpsi.EQ.O)  GOTO  82 

81  z=x2+rpsi*romega*(x2-xi)*exp(-romega/2/rpsi*tgap) 
GOTO  83 

82  z=0 

83  IF (z.GE. delta)  goto  20 
disp=2*x2-xl+c*f 2 
CALL  DATA 
tgap=tgap+dt 
GOTO  81 

100  CONTINUE 
END 

SUBROUTINE  DATA 

COMMON  J ,T,X1 .X2.F2, DT .ROMEGA . RPSI , DISP 
COMMON  RMASS, DELTA, BU.BL.CF 
j=j  +  l 
t=t+dt 

C  PRINTING  ALTERNATING  VALUES 
IF ( j . NE.2)  GOTO  400 
vel=x2-xl 

C  WRITE(unit=16,fmt=*)  t,X2 
WRITE (unit=17,fmt=*)  f2,X2 
j=0 

IF(t .GT. 131 . )  STOP 
C  GETTING  NEXT  DATA  POINT 
400  xl=x2 
x2=disp 

READ(unit=15 ,fmt=*)  dummy,  f2 
bu=delta-cf *(x2-xl) 
bl=-delta-cf *(x2-xl) 

RETURN 

500  CONTINUE 
END 
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C.2  The  RLS  Algorithm  ^ 

C  PROGRAM  "RLS" 

DIMENSION  U(3) ,  Y(3) .  A(6).HINV(6) 
DOUBLEPRECISION  X1.X2.X3.DEN 
OPEN (UNIT=17 , STATUS= *  OLD  * ) 

OPEN (UNIT=18 .STATUS3 ’NEW*) 

OPEN (UN IT=19 , STATUS= *  NEW ’ ) 

OPEN (UNIT=20 . STATUS3 • NEW ’ ) 

OPEN (UN IT=62 .STATUS3 ' OLD ’ ) 

K=0 

READ (UN IT=62 . FMT3*)  GAMMA , DT , NDATA 
SKIP=0 
A ( 1 ) =0 
A(2)=0 
A(3)=0 
A(4)=0 
A(5)=0 
A(6)=0 
FZ=0 
GZ=0 
HZ=0 
C 

READ (UNIT317, FMT3*)  U(1),Y(1) 

RE AD (UN I T= 1 7 , FMT3  * )  U(2),Y(2) 

DO  200  N=l, NDATA 

READ (UNIT3 17 , FMT3*)  U(3),Y(3) 

C 

Z=(Y (3) -2*Y (2) +Y(1)) /DT**2 
F=(Y(1)-Y(3))/DT/1000 
G*-Y (2) 

H3U(2) 

A(1)=A(1)+F**2 

A(2)=A(2)+F*G 

A(3)=A(3)+F*H 

A(4)=A(4)+G**2 

A(5)3A(6)+G*H 

A(6)=A(6)+H**2 


FZ*FZ+F*Z  *" 

HZ=HZ+H*Z 

GZ=GZ+G*Z 


C 

IF(K.LT.IO)  GOTO  190 

DEN=A(1) * (A (4) *A(6) -A(5) **2) +A(2) * (A (3) *A(5) -A (2) *A(6) ) 
DEN=DEN+A(3) *(A(2)*A(5)-A(3) *A(4) ) 

IFCDEN.EQ.O.)  GOTO  190 
HINV(1)= (A(4)*A(6) -A (5) * *2) /DEN 
HINV(2)-(A(3) *A(5)-A(2)*A(6))/DEN 
HINV(3)=(A(2)*A(6)-A(3)*A(4))/DEN 
HINV(4)=(A(1) *A(6) -A (3)* *2) /DEN 
HINV(6)=(A(2) *A(3)-A(i)*A(5) )/DEN 
HINV(6)=(A(1)*A(4)-A(2)**2)/DEN 
X1=HINV(1) *FZ+HINV (2) *GZ+HINV(3) *HZ 
X2=HINV(2) *FZ+HINV (4) *GZ+HINV (6) *HZ 
X3=HINV(3) *FZ+HINV(5) *GZ+HINV(6) *HZ 
C 

SKIP=SKIP+1 

IFCSKIP.LT. 6)  GOTO  190 
T=DT*N 


R0MEGA=ABS(X2)**.5 
RPSI=X1/R0MEGA/ 1000 
IFCX3.EQ.0.)  GOTO  190 
RMASS=1/X3 

WRITE (UNIT=18,FMT=*)  T.ROMEGA 
SI=RPSI*100 

WRITE(UNIT=19,FMT=*)  T,SI 
WRITE (UNIT=20,FMT=O  T.RMASS 
SKIP=0 


190  U(1)=U(2) 
U(2)=U(3) 
Y(1)=Y(2) 
Y(2)=Y(3) 
K=K+1 

200  CONTINUE 
END 


C.3  The  RLLS  Algorithm 


C  PROGRAM  "RLLS" 

C 

C  INITIALIZATION 

DOUBLEPRECISION  P(2,2.2) .PSTAR(2,2,2) ,S1(2.2) 
DOUBLEPRECISION  GMMAC2) , A,B , C , AS ,BS ,CS 
DIMENSION  H0LD(2 ,2) ,R(2,2) ,E(2,2) 

DIMENSION  Gl(2,2) ,H1(2,2) 

DIMENSION  G2 (2 . 2) , H2 (2 , 2) 

REAL  K(2,2.2) ,KSTAR(2.2,2) 

REAL  KBAR(2) ,KAST(2) 

DOUBLEPRECISION  DEN 
OPEN (UNIT-29 . STATUS=  *  OLD ’ ) 

OPEN (UNIT=25 .STATUS® ' OLD  * ) 

OPEN (UNIT=26 . STATUS=  *  NEW  * ) 

OPEN (UNIT=27 . STATUS® ' NEW  * ) 

OPEN (UNIT=62 . STATUS® ' OLD ’ ) 

C 

DO  30  12=1.2 
DO  20  1=1,2 
DO  10  >1,2 
K(I2,I, J)=0 
KSTAR(I2 , I , J) =0 
P(I2,I , J)=0 
PSTAR(I2,I, J)=0 
S1(I,J)=0 
10  CONTINUE 
H0LD(2,I)=O 
P (12. 1. 1)  =  100000 
PSTAR(I2 , I , I) =100000 
20  CONTINUE 
30  CONTINUE 
c 

READ (UN IT=52 , FMT=* )  GMMA(l) .DT.ndata 

AMBDA-1.0 

KSKIP=0 

C  NOTE:  GMMA® (L) AMBDA*VARIANCE 


I 


READ (UN IT=29 , FMT=  * )  HOLD (2, 2) .HOLD (2,1) 

READ (UN IT=29 , FMT* * )  R(1.2) ,R(1,1) 

DO  800  ITIME^l .NDATA 
READ(UNIT=29,FMT=*)  E(1.2) .E(l.l) 

C  RECURSION  0RDER=2 
DO  400  N=1 ,2 
C  FACTORIZE  P  TO  UDUt 
B=P(N,2,2) 

C=P(N ,2 , 1) /B 

A=P(N , 1 , 1)-P(N,2, 1)**2/B 
C  CALCULATE  kbar 

DEN~A*R(N . 1) **2+B* (C*R(N . 1) +R(N ,2) ) **2+GMMA(N) 
IF(DEN.NE.O)  GOTO  31 
DEN=. 00000001 

31  KBAR(l) = (P(N , 1 , 1) *R(N . 1) +P(N . 1 , 2) *R(N .2) ) /DEN 
KBAR(2)=(P(N,2.1)*R(N,1)+P(N,2,2)*R(N,2))/DEN 
C  CALCULATE  K(t) 

T1=E(N , 1) - (K(N , 1 , 1)*R(N , 1)+K(N , 1 ,2)*R(N ,2) ) 

T2=E(N ,2) - (K(N ,2 , 1) *R(N, 1) +K(N , 2 ,2) *R(N ,2) ) 

C  LINE  60 

K(N,1,1)=K(N,1, 1) +T1*KBAR(1) 

K(N ,  1 ,2)=K(N, 1 ,2)+Tl*KBAR(2) 

K(N .2 . 1)=K(N ,2 , 1) +T2*KBAR(1) 

K(N ,2,2)=K(N ,2 ,2) +T2*KBAR(2) 

C  ‘  CALCULATE  P(t-l) 

Z1=(P(N , 1 , 1) *R(N , 1) +P(N , 1 ,2) *R(N , 2) ) **2 
Z2=(P(N ,1,1)*R(N,1)+P(N,1,2)*R(N,2)) 

Z2=Z2* (P(N, 1 ,2) *R(N . 1) +P (N . 2 . 2) *R(N .2) ) 
Z4=(P(N.1,2)*R(N,1)+P(N.2,2)*R(N,2))**2 
P(N ,  1 , 1)  =  (P(N , 1 , 1) -Z1/DEN)/AMBDA 
P(N , 1 ,2)=(P(N , 1 ,2) -Z2/DEN)/AMBDA 
P(N ,2, 1)=P(N , 1 ,2) 

P(N ,2 ,2)=(P(N ,2 ,2) -Z4/DEN)/AMBDA 
C  CALCULATE  KAST 

BS=PSTAR(N ,2,2) 

CS=PSTAR(N,2, 1)/BS 

AS=PSTAR(N . 1 , 1) -PSTAR(N , 2 , 1) **2/BS 

DEN=AS*E(N . 1) **2+BS* (CS*E(N . 1) +E(N ,2) ) **2+GMMA(N) 


IF(DEN.NE.O)  GOTO  128 
DEN “.00000001 

128  KAST ( 1) “ (PSTAR(N .  1 , 1)  *E  (N  .  1) +PSTAR (N ,  1 , 2)  *E  (N ,  2)  ) /DEN 
KAST(2)“(PSTAR(N ,2,1) *E(N ,1) +PSTAR(N ,2 ,2) *E(N ,2))/DEN 
C  CALCULATE  KSTAR(t) 

T1=R(N , 1) - (KSTAR(N , 1 , 1) *E(N , 1) +KSTAR(N , 1 ,2) *E(N ,2) ) 

T2=R(N ,2) - (KSTAR(N , 2 , 1) *E(N , 1) +KSTARCN , 2 , 2) *E(N .2) ) 
KSTAR(N,1 , 1)=KSTAR(N,1 , 1)+T1*KAST(1) 

KSTAR(N ,1,2) =KSTAR(N ,1.2) +T1*KAST(2) 

KSTAR(N ,2 , 1) =KSTAR(N , 2,1) +T2*KAST(1) 

KSTAR (N .2.2) =KSTAR(N .2.2) +T2*KAST (2) 

C  CALCULATE  PSTAR(t) 

Yl= (PSTARCN , 1 , 1) *E(N . 1) +PSTARCN , 1 . 2) *E(N . 2) ) **2 
Y2= (PSTARCN . 1 . 1) *E(N , 1) +PSTAR(N . 1 , 2) *E(N . 2) ) 

Y2=Y2* (PSTARCN , 1 . 2) *E (N , 1) +PSTAR(N . 2 , 2) *E (N , 2) ) 

Y4= (PSTARCN . 1 . 2) *E (N , 1) +PSTAR(N ,2.2)*E(N,2))**2 
PSTARCN . 1 . 1)“ (PSTARCN . 1 . 1) -Yl/DEN) /AMBDA 
PSTARCN .1.2)= (PSTARCN .1.2) -Y2/DEN) /AMBDA 
PSTARCN .2 , 1)=PSTAR(N .1.2) 

PSTARCN .2.2)= (PSTARCN ,2,2) -Y4/DEN) /AMBDA 
IF(N.GT.l)  GOTO  400 
C  UPDATE  ORDER 
M=N*1 

DO  160  1=1,2 
DO  130  J-1,2 
G1  (I , J)=K(1 ,1, J) 

HI (I , J)=KSTAR(1 . I , J) 

130  CONTINUE 

E(M, I)=E(N . I) - (K(N , I , 1) *R(N , 1) +K (N , I ,2) *R(N ,2) ) 

R(M, I)=HOLD(M, I) 

HOLDCM, I)=R(N , I) ~(KSTAR(N , I , 1)*E(N , 1)+KSTAR(N. I ,2)*E(N  2)) 
150  CONTINUE 
C  LINE  100 

B=P(N,2,2) 

C=P(N,2,1)/B 

A=P(N . 1 , 1)-P(N,2, 1) **2/B 

T3=A*R(N .1)**2+B*(C*R(N,1)+RCN,2))**2 

GMMA(M)=GMMA(N) -AMBDA+T3 


C  CLOSE  LOOP  FOR  2ND  STAGE 
400  CONTINUE 

C  CALCULATE  PARAMETER  MATRICES 
DO  550  1-1,2 
DO  500  J=1 ,2 

G2(I ,  J)=G1  (I .  J) ~(K(2,I, 1)*S1(1 , J'  K(2 , I ,2) *S1 (2 , J) ) 

H2 ( I , J ) =KSTAR (2 , 1 ,  J  ) 

500  CONTINUE 

DO  510  J=3 ,4 
L=J-2 

G2(I, J)=K(2,I,L) 

H2(I , J)=S1 (I ,L) - (KSTAR(2 , I , 1) *G1 (1 ,L) +KSTAR(2 ,1,2) *G1 (2 ,L) ) 
610  CONTINUE 
550  CONTINUE 
C  RECONSTRUCT  PARAMETERS 
C  SOT=PSI*OMEGA*DT 

KSKIP  *  KSKIP+1 
IF (KSKIP.LT. 6)  GOTO  695 
S0T=(1+G2(1 ,3))/(l-G2(l ,3)) 

ROMEGA=ABS (2-G2 ( 1 , 1 ) * ( 1 +SOT) ) *  * . 5/DT 
RPSI=SOT/ROMEGA/DT 
IF(G2(1,2) .EQ.O)  GOTO  695 
IF(SOT.Eq.-l)  GOTO  695 
TME=DT*ITIME 

RMASS* (DT**2) / (G2 ( 1 , 2) +G2 ( 1 . 4) ) / ( 1+SOT) 

C  WRITE (UN IT=25,FMT=*)  TME.ROMEGA 

SI=RPSI*100 

C  WRITE(UNIT=26,FMT=*)  TME.SI 

C  WRITE (UNIT=27 ,FMT=*)  TME.RMASS 

KSKIP=0 

C  SHIFT  ANY  TIME  PARAMETERS 
695  DO  710  1*1 ,2 
DO  700  J=1 , 2 
SI(I.J)-HKI.J) 

700  CONTINUE 

R(1.I)*E(1.I) 

710  CONTINUE 
C  NEW  DATA  POINT 


800  CONTINUE 

WRITE (UNIT=25 . FMT=*)  ROMEGA . SI , RMASS 
END 


o  o 


C.4  The  Extended  RLS  Algorithm 

PROGRAM  "ERLS” 

DIMENSION  A(10).B(10).  U(3),  Y(3) 
OPEN (UNIT=17 , STATUS= *  OLD ' ) 

OPEN (UNIT=21 , STATUS3 ‘ NEW ’ ) 

OPEN (UNIT=22 , STATUS3 '  NEW  * ) 

OPEN (UN IT=23 , STATUS3  *  NEW ’ ) 

OPEN (UNIT=24 . STATUS3  *  NEW 1 ) 

C 

NDATA=1251 
DT=. 10 
BU= .  1 
BL=- .  1 
J=-l 
C 

DO  60  1=1,10 
A(I)=0 
60  CONTINUE 
EZ=0 
FZ=0 
GZ=0 
HZ=0 

q=o 

H=0 

1=0 

K=0 

C 

100  READ (UN IT= 17 , FMT=  * )  U(l) ,Y(1) 

IF(Y(1) .GT.BU)  GOTO  110 
IF(Y(1) .LT.BL)  GOTO  120 

Q=qn 

GOTO  100 
110  I=-l 

GOTO  160 
120  1=1 

160  READ(UNIT=17 ,FMT3*)  U(2) ,Y(2) 


IF(Y (2) .GT.BU)  GOTO  160  >' 

IF(Y(2) .LT.BL)  GOTO  170 
166  CONTINUE 

q=q+2 

1*0 

GOTO  100 

1 60  IF(I.NE.-l)  GOTO  165 
GOTO  200 

170  IF(I.NE.l)  GOTO  155 
200  CONTINUE 
M=Q+3 
N=0 

DO  600  N=M , NDATA 

READ (UNIT=17 , FMT=*)  U(3),Y(3) 

IF (Y (3) .GT.BU)  GOTO  210 
IF (Y (3) .LT.BL)  GOTO  220 
205  CONTINUE 
q=N 
1=0 

GOTO  100 

210  IF(I.NE.-l)  GOTO  205 
GOTO  230 

220  IF(I.NE.l)  GOTO  205 

230  Z=- (Y (3) -2*Y (2) +Y( 1) ) /DT**2 
E=(Y (3)-Y(l))/DT/3 
F=Y (2) 

G=-U(2)/2.0 
H=I*2 .0 

260  A(1)=A(1)+E**2 
A(2)=A(2)+E*F 
A(3)=A(3)+E*G 
A(4)=A(4)+E*H 
A(5)=A(6)+F**2 
A(6)=A(6)+F*G 
A(7)=A(7)+F*H 
A(8)=A(8)+G**2 


A(9)*A(9)+G*H 
A(10)«A(10)+H**2 
DO  260  L=1 , 10 
IF(A(L) .NE.O)GOTO  260 
A(L)*1.E-10 
260  CONTINUE 
C 

EZ=EZ+E*Z 

FZ=FZ+F*Z 

GZ=GZ+G*Z 

HZ=HZ+H*Z 

C  INVERTING  HTH  (4x4) 

C  USEFUL  VARIABLES: 

C1=A(1)*A(5)-A(2)**2 

C2=A(2)*A(6)-A(3)*A(6) 

C3=A(3)*A(7)-A(4)*A(6) 

C4=A(6)*A(9)-A(7)*A(8) 

C5=A(8)*A(10)-A(9)**2 

C6=A(6) *A(10)-A(9)*A(7) 

C 

D1=A(1) *A(6) -A (2) *A(3) 

D2=A(2) *A(7) -A (4) *A(5) 
D3=A(1)*A(7)-A(2)*A(4) 

D4=A(3)*A(9) -A(4)*A(8) 

D5=A(3) *A(10) -A (4) *A(9) 

C  COMPUTING  INVERSE 

IFCK.LT.4)  GOTO  230 

DEN=C1*C5-D1*C6+D3*C4+D5*C2-D4*D2+C3**2 
IF(DEN.EQ.O)  GOTO  280 
B(1)=(A(5)*C5-A(6) *C6+A(7) *C4)/DEN 
B(2)=(A(3)*C6-A(2) *C5-A(4)*C4)/DEN 
B(3)=(A(10) *C2-A(9) *D2+A(7)*C3)/DEN 
B(4)=(A(8)*D2-A(9)*C2-A(6)*C3)/DEN 
B(6)=(A(1) *C5+A(4)*D4-A(3) *D5)/DEN 
C 

B(6)=(A(2)*D5-A(1)*C6-A(4) *C3)/DEN 
B(7)=(A(1)*C4-A(2) *D4+A(3)*C3)/DEN 
B(8)=(A(10)*C1-A(7) *D3+A(4) *D2)/DEN 


V3 
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B(9)=(A(7)*D1-A(4)*C2-A(9)*C1)/DEN 
B(10)-(A(8)*C1+A(3)*C2-A(6)*D1)/DEN 
C  CALCULATE  VARIABLES 

X1=B(1)*EZ+B (2) *FZ+B (3) *GZ+8 (4) *HZ 
X2=B(2)*EZ+B(6)*FZ+B(6) *GZ+B (7) *HZ 
X3=B(3) *EZ+B(6) *FZ+B (8) *GZ+B (9) *HZ 
X4=B(4)*EZ+B(7) *FZ+B(9) *GZ+B(10) *HZ 


Xl=Xl/3 

X3=X3/2.0 

X4=X4*2 .0 

R0MEGA=ABS(X2)**.5 

RPSI=X1/R0MEGA 

RMASS=1/X3 

DELTA® (X4/X2) 

C  CALCULATE  FOR  NEW  GAP  SIZE 
CF=2*RPSI/R0MEGA 
IF  (K.LT.200)  GOTO  270 
VEL=(Y(2)-Y(1))/DT 
IFCVEL.GT.O.)  GOTO  265 
BU=1 . 1* (DELTA-CF+VEL) 
BL=-DELTA 
GOTO  270 
265  BU=DELTA 

BL=1 . 1* (0-DELTA-CF*VEL) 

C 

270  T=N*DT 
J=J+1 

IF(J.LT.IO)  GOTO  280 

WRITE (UNIT=21 ,FMT=*)  T.ROMEGA 

SI=RPSI*100 

WRITE (UNIT=22,FMT=*)  T.SI 
WRITE (UNIT=23,FMT=*)  T.RMASS 
WRITE (UNIT=24.FMT=*)  T. DELTA 
J=0 

280  K=K+1 
C 

490  Y(1)=Y (2) 
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